Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation 
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We present a further theoretical extension to the kinetic theory based formulation of the lattice Boltzmann 
method of Shan et al (2006). In addition to the higher order projection of the equilibrium distribution function 
and a sufficiently accurate Gauss-Hermite quadrature in the original formulation, a new regularization proce- 
dure is introduced in this paper. This procedure ensures a consistent order of accuracy control over the non- 
equilibrium contributions in the Galerkin sense. Using this formulation, we construct a specific lattice Boltz- 
mann model that accurately incorporates up to the third order hydrodynamic moments. Numerical evidences 
demonstrate that the extended model overcomes some major defects existed in the conventionally known lattice 
Boltzmann models, so that fluid flows at finite Knudsen number (Kn) can be more quantitatively simulated. Re- 
sults from force-driven Poiseuille flow simulations predict the Knudsen's minimum and the asymptotic behavior 
of flow flux at large Kn. 

PACS numbers: 47.11,+j, 51.10.+y, 47.45.Gx, 47.85.Np 



I. INTRODUCTION 

Understanding and simulating fluid flows possessing sub- 
stantial non-equilibrium effects pose a long standing chal- 
lenge to fundamental statistical physics as well as to many 
other science and engineering disciplines Illy]. Due to either 
rarefaction effects or small geometric scales, such flows are 
characterized by a finite Knudsen number, defined as the ratio 
between the particle mean free path, I, and the characteristic 
length, L, Kn = l/L. At sufficiently large Knudsen numbers, 
many of the continuum assumptions breakdown |3]. In par- 
ticular, the Navier-Stokes equation and the no-slip boundary 
condition become inadequate. 

Since the Boltzmann equation is valid for describing fluid 
flows at any Kn |4], the conventional approach for construct- 
ing extended hydrodynamic equations for higher Kn regimes 
has been through employing higher order Chapman-Enskog 
approximations resulting in, e.g., the Burnett and super Bur- 
nett equations. However, this approach encounters both theo- 
retical and practical difficulties O0I. Alternatively, attempts 
have been made to extend the Grad's 13 moment system [7] 
by including contributions of higher kinetic moments [8j]. One 
major difficulty has been the determination of the boundary 
condition for these moments because only the lowest few have 
clear physical meanings. In addition, due to the complexity in 
the resulting equations, application of this approach is so far 
limited to simple one-dimensional situations. Nevertheless, 
the moment based formulation offers an valuable insight into 
the basic fluid physics for high Kn flows. 

Over the past two decades, the lattice Boltzmann method 
(LBM) has developed into an efficient computational fluid dy- 
namic (CFD) tool 1 9]. Due to its kinetic nature, LBM intrin- 
sically possesses some essential microscopic physics ingre- 
dients and is well suited for handling more general bound- 
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ary conditions. Certain characteristic phenomena in micro- 
channel flows were predicted in LBM simulations at least 
qualitatively C3 EH [S E3- Q EH EE O E3 - In addition, 
by introducing a "stochastic virtual wall collision" process 
mimicking effects of free particle streaming in a long straight 
channel 1 14], analytically known asymptotic behavior at very 
large Kn were also produced. Nevertheless, being historically 
developed only to recover fluid physics at the Navier-Stokes 
level, the existing LBM schemes used in these studies pos- 
sess some well known inaccuracies and numerical artifacts. 
Therefore, strictly speaking the schemes are not applicable to 
high Kn flows other than for some rather limited situations. It 
is important to develop an LBM method capable of perform- 
ing accurate and quantitative simulations of high Kn flows in 
general. 

Recently, based on the moment expansion formulation 1 19], 
a systematic theoretical procedure for extending LBM beyond 
the Navier-Stokes hydrodynamics was developed 1 20]. In this 
work, we present a specific extended LBM model from this 
procedure containing the next order kinetic moments beyond 
the Navier-Stokes. A three-dimensional (3D) realization of 
this LBM model employs a 39-point Gauss-Hermite quadra- 
ture with a sixth order isotropy. In addition, a previously re- 
ported regularization procedure 121112211 . that is fully consis- 
tent with the moment expansion formulation, is incorporated 
and extended to the corresponding order. Simulations per- 
formed with the extended LBM have shown to capture certain 
characteristic features pertaining to finite Kn flows. There is 
no empirical models used in the new LBM. 



D. BASIC THEORETICAL DESCRIPTION 

It is theoretically convenient to describe a lattice Boltz- 
mann equation according to the Hermite expansion represen- 
tation 1 20]. The single-particle distribution functions at a set 
of particular discrete velocity values, {£ a : a = 1, • • • , d}, 
are used as the state variables to describe a fluid system. The 
velocity-space discretization is shown to be equivalent to pro- 
jecting the distribution function onto a sub-space spanned by 



the leading N Hermite ortho-normal basis, denoted by M N 
hereafter, provided that {£ a } are the abscissas of a sufficiently 
accurate Gauss-Hermite quadrature I19t 12011 . Adopting the 
BGK collision model |23], the discrete distribution values, f a , 
satisfy the following equation: 
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where r is a relaxation time, /„ ' is the truncated Hermite ex- 
pansion of the Maxwell-Boltzmann distribution evaluated at 
£ a , and F a is the contribution of the body force term. The 
truncation level determines the closeness of the above equa- 
tion to approximate the original continuum BGK equation. 
A Chapman-Enskog analysis reveals that the Navier-Stokes 
equation is recovered when the second order moment terms 
are retained. As the higher order terms are included, physical 
effects beyond the Navier-Stokes can be captured systemati- 
cally 0. 

In this work we use a specific model of Eq. JTal that consists 
of moments up to the third order, one order higher than the 
Navier-Stokes hydrodynamics in the conventional LBM mod- 
els 1 9]. For our present investigation of flows at high Kn but 
low Mach numbers (Ma), here we set the temperature, T, to 
a constant for simplicity. Denoting the local fluid density and 
velocity by p and u, and defining u a — £ a • u for brevity, in 
the dimensionless units in which all velocities are normalized 
by the sound speed (i.e., \/T - h ' 
compact form: 



1), f a takes the following 
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where u 2 = u ■ u, and w a is the quadrature weight corre- 
sponding to the abscissa £ a . The last term inside the brackets 
represents the contribution from the third-order kinetic mo- 
ments 1 24] which was shown to be related to the velocity- 
dependent viscosity [25] but generally neglected in the con- 
ventional lattice Boltzmann models. 



According to the previous analysis |20], the Gauss-Hermite 
quadrature employed for solving a third-order truncated sys- 
tem must be accurate with the polynomials up to the sixth or- 
der. For ease in implementing LBM model on Cartesian coor- 
dinates, we use the 3D Cartesian quadrature E^ 9 7 of Ref. 1 20]. 
Its 2D projection gives a quadrature E 2 \. The abscissas and 
weights of both quadratures are provided in Table U Both 
LBM models can be verified to admit isotropy for tensors of 
the form ^ w a£a ■ ■ ■ £a U P t° me sixth order instead of fourth 
in the conventional LBM models. Explicitly speaking, recov- 
ery of correct hydrodynamic physics up to the third order re- 



quires up to six-order isotropy conditions as follows |2* 
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where the roman subscripts i, j, . . ., n denote Cartesian com- 
ponents. In the above, &y is the Kronecker delta function, 

while S\, kl and S\A lrnn represent, respectively, the forth order 
and the sixth order generalizations: 
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Indeed, direct verification shows that these are satisfied by 
the 3D 39-speed model (S| 9 7 ) and its 2D projection (E 21 7 ). 
The conventional LBM schemes only satisfy the above mo- 
ment isotropy conditions up to the forth order (i.e., the Navier- 
Stokes order). It is also important to mention that there exist 
a few other lattice Boltzmann models satisfying sixth order 
isotropy 126, 12711. However they contain more discrete speeds 
and more complicated coefficients, and are difficult to extend 
to even higher orders. 

With the Cartesian quadrature above, Eq. ( fTal can be di- 
rectly discretized in physical space and time, yielding a stan- 
dard lattice Boltzmann equation of the form: 



fa(x + ta,t+l) = f a (x,t)-- f a {x,t) - f a ' 



p(0) 
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As usual, the "lattice convention" with unity time increment 
is used here. 



III. THE REGULARIZATION PROCEDURE 

An LBM computation is generally carried out in two steps: 
the streaming step in which f a at x is moved to x + £ a , and 
the collision step in which f a (x) is replaced with the right- 
hand-side of Eq. (|3- When viewed as a projection of the con- 
tinuum BGK equation into M. N , this dynamic process intro- 
duces an error due to the fact that f a does not automatically 
lie entirely within M N . Borrowing the language from spec- 
tral analysis, this is analogous to the aliasing effect. When the 
system is not far from equilibrium, such an error is small and 
ignorable. On the other hand, this error can be resolved via 
an extension of the "regularization procedure" previously de- 
signed for improvement in stability and isotropy I2lll22il . In 



TABLE I: Degree-7 Gauss-Hermite quadratures on Cartesian grid. 
Listed are the number of points in the symmetry group, p, abscissas, 
£ a , and the weights w a - Quadratures are named by the convention 
E% n where the superscript d and subscripts D and n are respec- 
tively the number of abscissas, dimension, and degree of algebraic 
precision. The subscript FS denotes permutations with full symme- 
try. Note that since all velocities are normalized with sound speed, 
the Cartesian grid spacing has a unit velocity of r — \/3/2. 
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terms of the Hermite expansion interpretation, the regulariza- 
tion procedure is more concisely described as the following. 
We split the post-streaming distribution into two parts: 
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where f' a is the deviation from the truncated Maxwellian, or 
the non-equilibrium part of the distribution. As f a already 
lies entirely in the subspace H , the projection is to ensure 
that the non-equilibrium contribution also lies in the same sub- 
space for all times, and only needs to be applied to f' a . Effec- 
tively, the projection serves as a filtering (or "de-aliasing") 
process to ensure the system stay inside the defined subspace 
M N in a Galerkin interpretation. 

The projection is to convert f' a to a new distribution f a 
which lies within the subspace spanned by the first three Her- 
mite polynomials. Using the orthogonality relation of the 
Hermite polynomials and the Gauss-Hermite quadrature, f a 
is given by the pair of relations: 



3 1 
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where 'HS n ' is the standard n-th Hermite polynomial 119112: 
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H^ k (i) = SiSjSk-SiSjk-SiSik-bSij, (8d) 

and a^ the corresponding Hermite expansion coefficient, 
both rank-n tensors. The first two Hermite coefficients vanish 
due to the vanishing contribution from the non-equilibrium 
distribution to mass and momentum. The second and third 
Hermite coefficients are: 
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where a^ 2 ' is traceless due to the conservation of energy. 
Clearly from the above construction or from direct verifica- 
tion, the projected distribution f' a gives the same second order 
(momentum) and third order fluxes as the original f a , 
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This is an essential step to preserve the required non- 
equilibrium properties affecting macroscopic physics. Fur- 
thermore, unlike f' a in which all higher order moments are in 
principle present, the projected distribution f a can be shown 
via the orthogonality argument to give zero contributions 
to fluxes higher than the defined third order above. Con- 
sequently, its physical implication is rather apparent: The 
regularization procedure filters out the higher order non- 
equilibrium moments that contain strong discrete artifacts due 
to the insufficient support of the lattice basis. 

Overall, given the discrete non-equilibrium distribution, its 
projection in H 3 is fully specified by: 
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Incorporating the regularization procedure, Eq. Q is modified 
to become 



fa(x + (ia,t+l)= /W + ( 1 - -) f a + F a 
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It is revealing to realize that the right-hand-side represents a 
damping operator acting on the "non-equilibrium" part of the 
dsitribution function. It is an immediate extension to assign 
a different relaxation time to each individual Hermite mode. 
Namely we can recast the collision operator into an equivalent 
yet slightly more general form: 
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FIG. 1 : (Color online) Peak velocity of a decaying shear waves as 
simulated by the 9-state (9s) and the 21-state (21s) standard BGK 
and regularized (REG) models. For each model, simulation is carried 
out with the wave vector aligned with either the lattice links or their 
diagonals. The latter is denoted by the post-fix "diag" 



where the linear matrix operator takes the following generic 
form 



7, ?b?6 + y 2 -p. < 



Mab = 5ab+[l- ~ 
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The additional coefficients, B\ and 82 can have values between 
zero and one separately. Indeed, when 8\ = 82 = 1, then 
we recover the result in ( II H . the third-order accurate projec- 
tion operator. On the other hand, comparing Eq. dl4i with 
the similar expression in Ref . 1 2 1 ] , we find that the latter is in 
fact a second-order projection operator (i.e., the first term in 
Eq. dl4t ). or simply 8\ = 1 and 62 = 0. It is to be empha- 
sized that the correct application of the third-order projection 
operator is based on the necessary condition of the third order 
Hermite quadrature. Hence it cannot be realized via conven- 
tional low order LBM such as the popular D3Q19 (or D2Q9). 

The explicit form of the body force term comes directly 
from the Hermite expansion of the continuum BGK equa- 
tion BCLc^l . Up to the third order, it can be expressed as: 
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To be noted here is that, whereas the first term is entirely due 
to the equilibrium part of the distribution, the term related 
to a/ 2 ) are contributions from the non-equilibrium part. To 
our knowledge, the non-equilibrium contribution in the body- 
force has not been explicitly considered in the existing LBM 
literature. Nevertheless, although it is expected to play an im- 
portant role at large Kn, at moderate Kn (< 1), no significant 
effects due to non-equilibrium contribution are observed in the 
numerical experiments in the present work. 



IV. NUMERICAL RESULTS 
A. Shear wave decay 

The first series of numerical simulations performed are for 
the benchmark problem of 2D sinusoidal shear wave decay. 
This set of tests is to evaluate impact from the increased or- 
der of accuracy and from the regularization procedure on the 
resulting isotropy. Obviously, a numerical artifact free wave 
decay should not depend on its relative orientation with re- 
spect to the underlying lattice orientation. The lattice orienta- 
tion dependent artifact has often plagued discrete fluid mod- 
els especially at finite Knudsen number when non-equilibrium 
effects are significant. For this purpose, we have defined 
two sets of simulations. In the first set, a sinusoidal wave 
with a wavelength (L) of 128 grid spacing is simulated on 
a 128 x 128 periodic domain. The initial velocity field is 
given by u x — uq sm(y/2TrL) and u y — 0. The wave vec- 
tor is aligned with the lattice. In the second set, the same 
sinusoidal wave is rotated by 45 degrees from the original ori- 
entation and simulated on a matching periodic domain size 
of 181 (= 128\/2) x 181. The Knudsen number, defined as 
Kn = 2tc s /L, is chosen to be 0.2, where r is the relaxation 
time and c s the sound speed. These two sets of simulations 
were conducted using four representative models: the popular 
2D 9-state (9-s) model (D3Q9), and the present 2D 21-state 
, £,d(21-s) model based on E%\, both with and without the reg- 
ulkrization process. Note, the 2D 9-state model only admits 
a second-order regularization projection as discussed in the 
preceding section. In discussions hereafter we shall refer the 
models without the regularization as the BGK models, and the 
ones with regularization the REG models. 

In Fig.n the dimensionless peak velocity magnitude, nor- 
malized by its initial value and measured at the 1/4 width of 
wavelength, is plotted against the non-dimensionalized time 
normalized by the characteristic decay time to = A 2 jv, where 
v = c 2 s (t — 1/2) is kinematic viscosity. As one can eas- 
ily notice, the decay rate of the shear wave for the 9-s BGK 
model is substantially different between the lattice oriented 
setup and the 45 degree one, exhibiting a strong dependence 
on the orientation of the wave vector with respect to the lat- 
tice. This indicates a strong anisotropy of the model at this 
Kn. This is consistent to our expectation that the 9-s BGK 
model was originally formulated to only recover the Navier- 
Stokes hydrodynamics {i.e., at vanishing Kn). Interestingly 
this anisotropy is essentially eliminated by the second-order 
regularization procedure in the resulting 9-s REG model. On 
the other hand, the amplitude of the shear wave exhibits a 
strong oscillatory behavior in addition to the exponential de- 
cay, implying a greater than physical ballistic effect. These 
may be explained as the following: The non-equilibrium part 
of the post-streaming distribution contains contributions from 
in principle all moments, which are highly anisotropic due 
to inadequate (only up to the second-order moment) sym- 
metry in the underlying discrete model. The regularization 
procedure filters out all the higher than the second order mo- 
ment contributions, yielding an isotropic behavior supported 
by the given lattice. On the other hand, the higher moments 



are critical at large Kn. Therefore, though isotropic, the 9-s 
REG model still should not be expected to show satisfactory 
physical results at high Knudsen numbers. For the 21-s BGK 
model, an anisotropic behavior is also very visible, though at 
a much smaller extent. This may be attributed to the residual 
anisotropy in the moments higher than the third order. Again, 
the anisotropy behavior is completely removed once the reg- 
ularization procedure is applied in the 21-s REG model, as 
shown from the totally overlapped curves between the lattice 
oriented and the diagonally oriented simulations. It is also 
noticeable that the decay history shows a much reduced oscil- 
latory behavior in the 21-s REG model. Because of its correct 
realization of the third order moment flux, we expect the re- 
sult is more accurate at this Knudsen number value. It is also 
curious to observe that the decay of the "lattice aligned" result 
from 9s-BGK is surprisingly close to that of 21-s REG model 
at this Kn. This is likely to be a mere coincidence. 



B. Finite Knudsen number channel flows 

Using the same four models, we subsequently carried out 
simulations of the force-driven Poiseuille flow for a range 
of Knudsen numbers. In order to identify the impact in ac- 
curacy in the resulting lattice Boltzmann models as oppose 
to the effects from various boundary conditions, here a stan- 
dard half-way bounce-back boundary condition is used in the 
cross-channel (y) direction. Furthermore, since we are not in- 
terested, for the present study, in any physical phenomenon 
associated with stream-wise variations, a periodic boundary 
condition is used with only four grid points in the stream-wise 
(x) direction. In the cross stream (y) direction, two different 
resolutions, L = 40 and 80 are both tried to ensure sufficient 
resolution independence. The Knudsen number is defined as 
Kn = v/(c s L). The flow is driven by a constant gravity force, 
g, pointing in the positive x direction. The magnitude of the 
force is set to SvUq/L 2 , which would give rise to a parabolic 
velocity profile with a peak velocity of Uq in the vanishing 
Kn limit. For consistence, a modified definition of fluid ve- 
locity, u — > u + g/2, is used in fa . Since the LBM models 
presently used here are all assuming constant temperature, to 
enforce an incompressible behavior with negligible thermo- 
dynamic effect throughout the simulated Kn range, we choose 
a sufficiently small value of Uq, corresponding to the nom- 
inal Mach numbers of Ma(= Uq/c s ) ~ 1.46 x 10~ 6 and 
2.92 x 10 -7 , and verified that our results are independent of 
Ma. The actual resulting fluid velocity in these simulations 
can achieve values much higher than Uq at higher Kn. 

Plotted in Fig. [2] is the non-dimensionalized mass flux, 
Q = J2„_ u x (y)/Qo, as a function of Kn in the final steady 
state of the simulations. Here Qo — gL 2 /c s . For compar- 
ison we also include two analytical asymptotic solutions [4] 
for both small and very large Kn. To be noted first is that at 
the vanishing Kn limit, all simulation results agree with each 
other as well as with the analytical solution. This confirms 
that all these LBM schemes recover the correct hydrodynamic 
behavior at vanishing Kn, i.e., the Navier-Stokes limit. Also 
plotted is the exact Navier-Stokes solution of Q = l/(12Kn), 



a well understood monotonically decreasing line with no min- 
imum. At higher Kn, the 9-s BGK model exhibits a Knudsen's 
minimum while over-estimates the flux according to some 
previously published reports 11411 . However, by filtering out 
moment contributions higher than the second order, such a 
phenomenon is completely disappeared from the result of 9- 
s REG, yielding a purely monotonically decreasing behavior. 
This is a rather interesting but not entirely surprising result. 
Once again, the regularization process enforces the system to 
be confined within the second-order Hermite moment space, 
while all higher order non-equilibrium contributions including 
both the numerical artifacts and the physical ones, responsible 
for the finite Knudsen phenomena, are filtered out. Conse- 
quently, only the Navier-Stokes order effects are preserved in 
the 9-s REG model no matter the degree of non-equilibrium 
at finite Kn. To be further noticed is the impact of the second- 
order regularization on the near wall properties. In particu- 
lar, Fig. [5] shows that 9-s REG gives vanishing slip velocity 
at the wall for a range of Kn values (Kn = 0.1, 0.2, and 
0.5). This suggests that a bounce-back boundary process is 
sufficient to realize a no-slip condition once a Navier-Stokes 
order dynamics in the model equation is enforced. This is 
yet another confirmation as to why the resulting Q from 9-s 
REG model lies very close to the exact Navier-Stokes theo- 
retical curve up to significantly high Kn values. In compar- 
ison, all the other three LBM models show finite slip veloc- 
ity values, due to the previous discussed reason that they all 
contain higher than the second-order non-equilibrium contri- 
butions. Once again, one must remember that the higher order 
effects in 9-s BGK (and to a lesser extent the 21-s BGK) have 
substantial lattice discrete artifacts. To be emphasized is the 
21-s REG model: Due to the regularization procedure, only 
the third-order non-equilibrium moment contribution is pre- 
served. On the other hand, because it has numerically shown 
to capture the Knudsen minimum phenomenon, correct incor- 
poration of the third-order moment physics is thus essential 
for accurately simulating some key flow physics beyond the 
Navier-Stokes. We also wish to emphasize that all these dif- 
ferences are due to the intrinsic natures in these LBM models, 
and has nothing to do with spatial and temporal resolutions. 
As a comparison, we plot in Fig. |4] velocity profiles generated 
from the 21-s REG model. 

There are a number of ways that gravity force can be in- 
cluded in LB equations. One can either treat the gravity force 
outside the collision operator as an external body force term as 
given by Eq. (1151 or by applying a local momentum/velocity 
shift inside the equilibrium distribution [30]. With regular 
BGK model, for all Kn, no difference in results are observed 
when the gravity force is applied in different ways. With 
the regularization procedure and sufficient isotropy however, 
some differences are observed at finite (> 1) Kn. Generally 
speaking, applying the body force via Eq. dl5> tends to pre- 
dict higher flux then via momentum/velocity shifting. For the 
results shown in Fig. [2] the velocity shift applied in fa is 
l/2g. The other half of the gravity force is applied as an ex- 
ternal body force term JBJ (Cf., J2J] ). 

The results from both the 21-s BGK and the 21-s REG mod- 
els predict a Knudsen minimum which resembles that of the 
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FIG. 3: (Color online) This figure shows the normalized velocity 
profiles at resolution 40 and Kn — 0.1, 0.2, and 0.5. Notice the 9-s 
REG model exhibits no visible slip velocity for both small and large 
Kn. 



9-s BGK except with reduced over-estimations at higher Kn. 
What is interesting, and requires further understanding, is that 
the flux behavior predicted by the 2 1 -s REG exhibits a reversal 
of curvature at higher Kn, resembling the analytical asymp- 
totic solution of Cercignani |4]. 

The qualitative differences seen from these four models 
suggest that contributions from moments beyond second or- 
der are essential for capturing fundamental physical effects 
at high-Kn. Although the high-order moments are implicitly 
contained in the second-order BGK model, its dynamics is 
highly contaminated with numerical lattice artifacts. In con- 
trast, by incorporating the high-order moments explicitly and 
systematically with the regularization, flows at these Kn val- 
ues can indeed be captured by the extended LBM model. 



FIG. 4: (Color online) This figure shows the normalized velocity 
profiles at resolution 40 and Kn = 0.1, 0.2, and 0.5. Notice the 21-s 
REG model exhibits substantial slip velocity at large Kn. 



V. CONCLUSION 

In summary, the kinetic based representation offers a well- 
posed approach in formulation of computational models for 
performing efficient and quantitative numerical simulations of 
fluid flows at finite Kn. In this work we present a specific 
extended LBM model that accurately incorporates the physi- 
cal contributions of kinetic moments up to the third order in 
Hermite expansion space. The new regularization procedure 
presented in this paper ensures that both the equilibrium and 
the non-equilibrium effects are confined in the accurately sup- 
ported truncated subspace at all times so that the unphysical 
artifacts are filtered out in the dynamic process. This resulting 
LBM model is robust and highly efficient. Because of its ac- 
curate inclusion of the essential third order contributions, this 
model is demonstrated to be capable of quantitatively captur- 
ing certain fundamental flow physics properties at finite Knud- 
sen numbers. This is accomplished without imposing any em- 
pirical models. Furthermore, because of the removal of dis- 
crete anisotropy, it is also clear that the new LBM model is 
not limited to specific uni-directional channel flows, nor it is 
only applicable for lattice aligned orientations. 

Nevertheless, a number of issues are awaiting further stud- 
ies. For even higher Kn (~ 10), one should expect moment 
contributions higher than the third order to become physically 
important. This is straightforward to include via the system- 
atic formulation 1 20] together with the regularization proce- 
dure described here. The issue of boundary condition is also 
of crucial importance 13 ll l32l 13311 . Even though, as demon- 
strated in this paper that the realization of the essential slip- 
velocity effect and the asymptotic behavior is attributed to a 
significant extent to the third order and higher moment con- 
tributions in the intrinsic LBM dynamic model itself. As re- 
ported in some previous works I13l UM . the kinetic bound- 
ary condition of Ansumali and Karlin 1 34] has led to substan- 
tial improvements at the Navier-Stokes level micro-channel 
flow simulations. Boundary condition itself serves as an ef- 



fective collision so that it also modifies the degree of the 
non-equilibrium distribution. Specifically, the well known 
Maxwell boundary condition (in which the particle distribu- 
tion function is assumed to be at equilibrium) should be ex- 
pected give different effect at high Kn compared with that of 
the bounce-back used in the present study. The latter bound- 
ary condition preserves non-equilibrium contributions at all 
orders. As a consequence, we suspect that the finite slip phe- 
nomenon is likely to still be over predicted in our current sim- 
ulations with the bounce-back boundary condition than what 
actually would occur in reality. Many more detailed and fur- 
ther investigations particularly pertaining channel flows at fi- 
nite Knudsen number are certainly extremely important in the 
future studies. 

Thermodynamic effect is also expected to become impor- 
tant at finite Kn when Mach number is not negligibly small. 
Some distinctive phenomena that is substantial and character- 
istic only at finite Knudsen numbers with sufficiently large 
Mach numbers 1 35 ] . The work of Xu 13511 demonstrates 



the importance of an accurate formulation of higher than 
the Navier-Stokes order thermodynamic effect at finite Kn. 
Theoretically, this additional property is associated with the 
so called super-Burnett effect in the more conventional lan- 
guage 1 35], or the forth-order moment contribution in the Her- 
mite expansion and can be incorporated in a further extended 
LBM model [20]. The present third-order model is thermody- 
namically consistent, but only at the Navier-Stokes level 1 20]. 
Another interesting point to mention is that both the 2 1 -s BGK 
or the 21-s REG can allow an expanded equilibrium distribu- 
tion form including terms up to the fourth power in fluid ve- 
locity (as opposed to the square power in 9-s BGK), for that 
the correct equilibrium energy flux tensor is still preserved. 
Including the forth power terms immediately results in a de- 
sirable positive-definite distribution for the zero-velocity state 
at all Mach number values. 
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